Identification of hub genes and drug candidates for NF2-related vestibular schwannoma by bioinformatics tools

Neurofibromatosis type 2 (NF2)-related vestibular schwannoma (NF2-VS) is a rare genetic disorder that results in bilateral acoustic neuromas. However, the exact pathogenesis of the disease is still unclear. This study aims to use bioinformatics analyses to identify potential hub genes and therapeutic. We retrieved the mRNA expression profiles (GSE108524 and GSE141801) of NF2-VS from the database, and selected the leading 25% genes with the most variance across samples for weighted correlation network analysis. Subsequently, we conducted gene ontology term and Kyoto Encyclopedia of Genes and Genomes signaling network enrichment analyses. The STRING database was employed for protein-protein interaction (PPI) axis construction. The mRNA-miRNA modulatory network was generated via the miRTarBase database. Differentially expressed genes (DEGs) were identified via the R package “limma” in both datasets, and hub genes were screened via intersection of common DEGs, candidate hub genes from the PPI axis, and candidate hub genes from the key module. Finally, common DEGs were uploaded onto the connectivity map database to determine drug candidates. Based on our observations, the blue module exhibited the most significant relation to NF2-VS, and it included the NF2 gene. Using enrichment analysis, we demonstrated that the blue modules were intricately linked to modulations of cell proliferation, migration, adhesion, junction, and actin skeleton. Overall, 356 common DEGs were screened in both datasets, and 33 genes carrying a degree > 15 were chosen as candidate hub genes in the PPI axis. Subsequently, 4 genes, namely, GLUL, CAV1, MYH11, and CCND1 were recognized as real hub genes. In addition, 10 drugs with enrichment scores < −0.7 were identified as drug candidates. Our conclusions offered a novel insight into the potential underlying mechanisms behind NF2-VS. These findings may facilitate the identification of novel therapeutic targets in the future.


Introduction
Vestibular schwannomas (VSs) are benign neoplasms of Schwann cell-derived vestibulocochlear nerves. [1,2]pproximately 90% of VSs patients experience progressive hearing loss, along with other symptoms like tinnitus and dizziness.The lifetime VS risk is estimated to be at 1 in 1000. [3]Most VSs appear sporadically, and are located primarily within the cerebellopontine angle.Approximately 5% of VS cases involve the neurofibromatosis type 2 (NF2), present bilaterally, and have a genetic predisposition. [4]Bilateral VS is a typical Nf2-related vestibular schwannoma (NF2-VS) lesion.The NF2-VS is a rare genetic disease brought on by the inactivation of the NF2 gene on chromosome 22q12.2. [4,5]The NF2 gene product, merlin, is generally expressed at high levels in Schwann and nerve cells among adult individuals.Given that merlin contributes heavily to cell-to-cell proliferation and adhesion, merlin deficiency can potentially impair contact growth inhibition function. [6]Till date, there are no reports on the exact mechanism of NF2-VS pathogenesis.
Gene expression profiling is a robust approach for identifying essential pathological mechanisms, [7] and one form of gene profiling, microarray-based gene expression profiling, is frequently used to screen for differentially expressed genes (DEGs)  or novel biomarkers in numerous human diseases. [7]The Gene Expression Omnibus database contains certain VS gene expression profiles, [8] the integration and reanalysis of which may provide valuable insight into the pathogenetic pathways of NF2-VS.Weighted correlation network analysis (WGCNA) is a bioinformatics tool that describes association profiles between different genes. [9]WGCNA incorporates genetic and clinical information to precisely identify relevance of certain genes among various samples.Herein, we, for the first time, screened for candidate NF2-VS-related hub genes, and identified several efficacious therapeutic drugs for NF2-VS using bioinformatics analysis.

Construction of a co-expression network
"WGCNA," [9] an R package, was used to conduct WGCNA.The "pickSoftThreshold" tool was employed for the softthresholding power measurement.Subsequently, a co-expression network module was developed using the "one-step network construction and module detection method," and the minimum gene quantity within the modules was set to 30.

Selection of the NF2-VS-related module
For the determination of major components of various modules, we employed module eigengenes (MEs), which utilize the first principal component of module expressions.We also evaluated the association between MEs and patient clinicophysiological information, and generated heat maps.

Functional enrichment analysis (FEA)
To elucidate the roles of NF2-VS-related genes within the modules, we performed FEA of the strongly associated modules using the clusterProfiler package of the R software. [10]P value adjustment was done via the Holm-Bonferroni technique.The adjusted P and the q value thresholds were set to 0.05.

Generation of a protein-protein interaction (PPI) axis
The key module genes were uploaded to STRING (https:// string-db.org/)to generate the PPI axis, and the results were visualized using the Cytoscape (v3.8.2) software. [11]Individual node degree was computed using Cytoscape.Lastly, genes with > 15 degree were considered to be candidate hub genes.

Construction and analysis of the miRNA-mRNA network
To explore correlations between miRNA and mRNA in the key module, we retrieved miRNA information from the miRTar-Base database, including, over 50,000 experimentally validated miRNA-target associations.The miRNA-mRNA modulatory network was visualized using Cytoscape.

Identification of drug candidates
To identify small molecular drugs that may effectively inhibit NF2-VSs occurrence and progression, we uploaded our list of up-and down-regulated genes into the connectivity map (CMap) database.CMap, a gene expression database, was employed to screen for functional associations among small molecular compounds, genes, and disease status. [13]A correlation score was then obtained from the CMap database based on the DEG enrichment in the reference gene expression.

Weighted gene co-expression network (WGCN) construction
The leading 25% of variant genes were entered in co-expression analysis.Sample hierarchical clustering methods were used to assess sample outliers in the GSE108524 and GSE141801 datasets, and GSM4213477 was eliminated from co-expression analysis (Fig. 1A and B).Next, the "pickSoftThreshold" tool from the WGCNA package was employed for soft-thresholding power (β) determination.The estimated GSE108524 power (β) was 7, however, no suitable power (β) was detected within the GSE141801 dataset (Fig. 1C and D).Thus, the WGCN was generated using 6305 genes from the GSE108524 dataset.Using WGCNA, 19 co-expression modules were identified (Fig. 1E  and F), and the gray modules represented genes that were not classified into other modules.

Identification of module-trait associations and key modules screening
MEs were assessed for individual modules, and the association between MEs and NF2-VSs was computed.The associated heatmap is displayed in Figure 2A.We observed a strong negative correlation between the blue module and NF2-VSs (R = −0.98 and P < .001),as depicted in Figure 2B.Additionally, the NF2 gene was also located within the blue module (Fig. 2A).Hence, the blue module was selected as the key module.Based on our WGCNA, 215 candidate hub genes were finally screened from the blue module.

FEA of genes in the blue module
Gene ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) FEA were conducted with the clus-terProfiler package, and our conclusions are presented in Figure 3.The analyzed genes were strongly enriched in the following biological process terms, namely, cell proliferation, cell migration, and movement; as well as the molecular function terms "mRNA binding involved in posttranscriptional gene silencing" and "RNA binding involved in posttranscriptional gene silencing," and cellular components related to the collagen-based extracellular matrix, cell-substrate junction, and focal adhesion (Fig. 3A).Moreover, the genes showed marked enrichment in 3 KEGG axes, namely, complement and coagulation cascades, tight junction, and actin cytoskeleton modulation (Fig. 3B).www.md-journal.com

PPI axis construction
The STRING database was used to generate the PPI axis of genes in the blue module.PPIs with interaction scores > 0.4 was employed for PPI axis construction.In all, 33 genes (15.35% of the candidate hub genes) with degrees > 15 were recognized as candidate hub genes belonging to the PPI axis (Fig. 4A).

The mRNA-miRNA modulatory network construction
The blue module consisted of 90 microRNAs (41.86% of the candidate hub genes).We retrieved 11 microRNAs (12.22% of the microRNAs) with |MM| >0.8 as candidate hub microRNAs.Subsequently, we employed the miRTarBase database to predict the target genes of 11 miRNAs.In all, 1586 genes were estimated, which were matched to 60 mRNAs (27.91% of the candidate hub genes) from the blue module.The miRNA-mRNA modulatory axis is displayed in Figure 4B.

Identification of DEGs and real hub genes
Following DEG screening based on the aforementioned criteria, 1124 and 1308 genes were identified from the GSE108524 and GSE141801 datasets, respectively.The top 100 DEGs for both datasets were showed in Figure 5A and C, respectively.These genes shared 356 DEGs (17.15% of total DEGs) between the 2 datasets, among which, 232 (65.17% of the shared DEGs) were scarcely expressed and 124 (34.83% of the shared DEGs) were highly expressed.The volcano plots for the DEGs of GSE108524 and GSE141801 datasets were presented in Figure 5B and D. Finally, 4 genes (1.12% of the shared DEGs) were screened as real hub genes upon intersection of the common DEGs, candidate hub genes from the PPI axis, and potential hub genes from the blue module (Fig. 5E and F).

Screening of potential therapeutic drugs
To explore drugs that potentially inhibit NF2-VSs development and progression, we converted the format of common elevated-and reduced gene expressions into probe IDs within the HG-U133A platform, and uploaded the information into the CMap database.Subsequently, negative enrichment scorebased stratification was used to identify leading small molecular compounds for drug development.Drug candidates with enrichment scores < −0.7 are provided in Table 1.

Discussion
NF2 is an autosomal dominant inherited disease brought on by variants within the NF2 gene.The classical NF2 symptom is bilateral VSs [14] ; the NF2 incidence at birth is approximately 1/33,000; and its overall prevalence is about 1/56,000. [15]NF2, harboring 17 exons, is a tumor suppressor gene.It encodes a protein named merlin containing 559 amino acids. [16]Merlin, belonging to the Ezrin/Radixin/Moesin family of membranecytoskeleton-interacting proteins, is primarily present within the plasma membrane and cytoskeletal compartments, and it interacts with multiple transmembrane receptors and intracellular proteins, namely, protocadherin fat (FAT), CD44, β1-integrin, and EGFR, thereby providing anchorage between the aforementioned proteins. [17]Despite strong evidence that merlin contributes to the stability of the membrane-cytoskeleton interface by suppressing the PI3kinase/Akt, Raf/MEK/ERK, and mTOR networks, [18][19][20] not much is known about the NF2 tumorigenesis process. [21]At present, bioinformatics analyses are used to elucidate the underlying mechanisms and potential therapeutic agents of numerous tumors.Herein, using WGCNA, we generated the NF2-based gene co-expression axes, and identified a key gene co-expression module strongly associated with NF2-VS pathogenesis.We also identified several drugs with potential therapeutic benefits for NF2-VSs using the CMap database.Based on our literature review, this study was the first to employ bioinformatics analysis to screen for hub genes and drug candidates for NF2-VSs.
Herein, we identified 19 co-expression modules using WGCNA.NF2-VS was most significantly correlated with the blue module.Moreover, since the NF2 gene was detected within the blue module, this module was then selected as the primary module.
GO and KEGG FEA are commonly used to elucidate the physiological roles of relevant genes.Herein, we employed the GO and KEGG analyses to examine the physiological activities of genes within the blue module.Based on the GO analysis, significant enrichments were found in cell proliferation and migration (in biological process), as well as cell adhesion and junction (in cellular components).Based on the KEGG network analyses, the actin cytoskeleton modulation was intricately related to NF2, which corroborated with the conclusions of prior studies.
In our attempt to identify real hub genes, we identified 356 common DEGs between NF2-VSs and normal nerves from the GSE108524 and GSE141801 datasets.We also employed STRING to construct the PPI axis of genes from the blue module.Ultimately, we identified 4 genes (GLUL, CAV1, MYH11, CCND1) as real hub genes via the intersection of common DEGs, genes with degree > 15 from the PPI axis, and potential hub genes from the blue module.GLUL encodes glutamine synthetase, the only enzyme that catalyzes the novel production of glutamine via condensation of ammonium and glutamate into glutamine. [22][25] Moreover, GLUL knockdown is reported to abrogate membranal localization and impair endothelial cell motility within human umbilical vein endothelial cells. [26]CAV1 is a membranal protein that coats the plasma membrane caveolae. [27]CAV1 phosphorylation is strongly associated with cellular transformation. [28,29]Prior investigations revealed that CAV1 reduces xenograft tumor formation from colon cancer cells within nude mice. [30][34][35] MYH encodes myosin heavy chain 11, a smooth muscle myosin that belongs to the myosin heavy chain family. [36]In recent years, a series of reports indicated that MYH11 modulates cell migration, association with cell adhesion proteins, and tumor suppression. [37,38]Moreover, mutations in the MYH11 gene is related to multiple cancer types, namely, colorectal, bladder, [39] and head and neck cancers. [40]Herein, CCND1 was the only hub gene that was elevated in NF2-VS samples.CCND1 encodes cyclin D1, a major cell cycle regulator. [41]Cyclin proteins D1 control cell cycle progression via modulating the G1-S phase transition. [42]CCND1 is an oncogene responsible for uncontrolled cell proliferation. [43]More importantly, merlin is known to inhibit CCND1. [44]ased on earlier reports, miRNAs are critical regulators of numerous tumors.Therefore, we selected miRNAs with high |MM| to generate a potential miRNA-targeted modulatory network.According to our analysis, miR-30c-1 (degree = 19) and miR-145 (degree = 13) are potential modulators of NF2-VSs.
We also identified candidate therapeutic drugs that may target common highly-and scarcely-regulated genes in NF2-VSs using the CMap database (Table 1).Proscillaridin, [45] naftopidil, [46] fursultiamine, and ursolic acid [47] are reported to possess robust antitumor activity.Among them, proscillaridin yields a higher enrichment score (enrichment = −0.928,P = .00058).Proscillaridin, a cardiac glycoside, is derived from Scilla and in Drimia maritima plants, and it possesses potent cytotoxic and anticancer properties.Collectively, our analyses indicated that these candidate drugs may be effectively employed to manage NF2-VSs.

Conclusions
In conclusion, herein, demonstrated that the key gene coexpression module, hub genes, and several functional networks, such as, actin cytoskeleton regulation, cell proliferation, and migration were strongly associated with NF2 pathogenesis.We also identified potential drugs that may inhibit NF2-VSs.These results are crucial for the future development and management of NF2-VSs.We warrant additional investigations into the underlying mechanism of hub genes and functional networks that potentially contributes to NF2 development.

Figure 1 .
Figure 1.Construction of a co-expression network.(A) Sample clustering dendrogram of GSE141801 using WGCNA.(B) Sample clustering dendrogram of GSE108524 using WGCNA.(C) The scale-free fit index assessment for several soft-threshold powers from the GSE108524 dataset.Red line adjusted to 0.85.(D) The scale-free fit index assessment for several soft-threshold powers from the GSE141801 dataset.Red line adjusted to 0.85.(E) The cluster dendrogram of the co-expression network-based genes.(F) The gene quantities in different modules.WGCNA = weighted correlation network analysis.

Figure 2 .
Figure 2. Screening of key modules.(A) Module-trait associations within the generated network.The top figure in each row denotes the relation to various clinical features, whereas, the bottom figure denotes the P value.The red box represents the module with NF2.(B) The NF2 gene relevance in the blue module (one dot refers to a single gene within the blue module).NF2 = neurofibromatosis type 2.

Figure 3 .
Figure 3. GO and KEGG network enrichment analyses of key module genes.(A) GO analysis of blue module genes.Only the top 20 BP terms are displayed.(B) KEGG network analysis of blue module genes.BP = biological process, GO = gene ontology, KEGG = Kyoto Encyclopedia of Genes and Genomes.

Figure 4 .
Figure 4.The blue module network.(A) PPI network of blue module genes, the node size indicates the network degree, the log 2 FC rang is presented as a color bar located at the bottom of the figure.The outer layer nodes represent genes with degrees > 15. (B) microRNA-target modulatory axis for the blue module.Blue triangles denote microRNAs, and green nodes denote genes.PPI = protein-protein interaction.

Figure 5 .
Figure 5. Identification of real hub genes.(A) Heatmap of top 100 DEGs from the GSE108524 dataset.(B) Volcano plot of DEGs in the GSE108524 dataset.(C) Heatmap of top 100 DEGs from the GSE141801 dataset.(D) Volcano plot of DEGs in the GSE141801 dataset.(E) Three genes were identified by intersecting the scarcely-expressed common DEGs, candidate hub genes from the PPI axis, and potential hub genes from the blue module.(F) One gene from the intersection of highly-expressed common DEGs, candidate hub genes from the PPI axis, and potential hub genes from the blue module.DEGs = differentially expressed genes, PPI = protein-protein interaction.

Table 1
Drug candidates for NF2.